home *** CD-ROM | disk | FTP | other *** search
/ Tech Arsenal 1 / Tech Arsenal (Arsenal Computer).ISO / tek-04 / bipl.zip / PROCS.ZIP / POM.ICN < prev    next >
Text File  |  1992-11-20  |  2KB  |  59 lines
  1. ############################################################################
  2. #
  3. #    File:     pom.icn
  4. #
  5. #    Subject:  Procedure to compute phase of the moon
  6. #
  7. #    Author:   Ralph E. Griswold
  8. #
  9. #    Date:     September 6, 1992
  10. #
  11. ###########################################################################
  12. #
  13. #  pom(n, phase) returns record with the Julian Day number of fractional
  14. #  part of the day for which the nth such phase since January, 1900.  Phases
  15. #  are encodes as: 0 - new moon, 1 - first quarter, 2 - full moon, 3 - last
  16. #  quarter.  GMT is assumed.
  17. #
  18. ############################################################################
  19. #
  20. #  Acknowledgement:  This procedure is based on an algorithm given in
  21. #  "Numerical Recipes; The Art of Scientific Computing"; William H. Press,
  22. #  Brian P. Flannery, Saul A. Teukolsky. and William T. Vetterling;
  23. #  Cambridge University Press, 1986.
  24. #
  25. ############################################################################
  26.  
  27. record jdate(number, fraction)
  28.  
  29. procedure pom(n, nph)
  30.    local i, jd, fraction, radians
  31.    local am, as, c, t, t2, extra
  32.  
  33.    radians := &pi / 180
  34.  
  35.    c := n + nph / 4.0
  36.    t := c / 1236.85
  37.    t2 := t * t
  38.    as := 359.2242 + 29.105356 * c
  39.    am := 306.0253 + 385.816918 * c + 0.010730 * t2
  40.    jd := 2415020 + 28 * n + 7 * nph
  41.    extra := 0.75933 + 1.53058868 * c + ((1.178e-4) - (1.55e-7) * t) * t2
  42.  
  43.    if nph = (0 | 2) then
  44.       extra +:=  (0.1734 - 3.93e-4 * t) * sin(radians * as) - 0.4068 *
  45.           sin(radians * am)
  46.    else if nph = (1 | 3) then
  47.       extra +:= (0.1721 - 4.0e-4 * t) * sin(radians * as) - 0.6280 *
  48.         sin(radians * am)
  49.    else fail
  50.  
  51.    if extra >= 0 then i := integer(extra)
  52.    else i := integer(extra - 1.0)
  53.    jd  +:=  i
  54.    fraction := extra - i
  55.  
  56.    return jdate(integer(jd), fraction)
  57.  
  58. end
  59.